#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Unbend I                                                           #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 02/20/2006                                             #
# Last Modification: 02/20/2006                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#  
# SORTORDER: 40
#
# MANUAL: The UNBEND I script performs the first unbending cycle. 
#
# MANUAL: A reference is generated by Fourier-filtering with a very small radius (holea, 1 pixel) and then windowing a small reference (boxa1, the width should be 10 times smaller than the image width). In parallel the original image FFT is masked with a larger radius (maska, ~20 pixel, 5 times smaller than the lattice spacing in Fourier space). This masked FFT and the Fourier-transformed and complex-conjugated reference are then multiplied (corresponding to a cross-correlation between image and reference), and the resulting cross-correlation map is searched for peaks (MRC program QUADSERCH). The MRC program CCUNBEND then unbends the original image, and the resulting unbent image is Fourier transformed to give an FFT with sharper diffraction peaks.
#
# MANUAL: The unbent FFT is evaluated (MRC program MMBOX), yielding a list of amplitudes and phases for each lattice spot, together with an IQ value for each spot. The statistic of IQ spots and the final average FFT peak profile is combined into one single value QVal, which can be used to judge the quality of the results of image processing.
#
# MANUAL: In this case, the actual script is not in the usual script directory, but in a separate script ${proc_2dx}2dx_unbend1_sub.com, which is called via the "source" command from the c-shell. For this reason, double-clicking the script name in the "Standard Scripts" window will only show the less interesting part of the script. The remainder can be found by using the editor directly to open the 2dx_unbend1_sub.com file in the 2dx_image/proc directory. This is designed such that it is easy to call the same core of the unbend1 job from the refinement scripts.
#
#
# PUBLICATION: 2dx - User-friendly image processing for 2D crystals: <A HREF="http://dx.doi.org/10.1016/j.jsb.2006.07.020">J. Struct. Biol. 157 (2007) 64-72</A>
# PUBLICATION: Image processing of 2D crystal images: <A HREF="http://link.springer.com/protocol/10.1007%2F978-1-62703-176-9_10">Methods in Molecular Biology (2012)</A>
#
# DISPLAY: maska
# DISPLAY: boxa1
# DISPLAY: boxa2	
# DISPLAY: quadrada
# DISPLAY: facthresha
# DISPLAY: quadpreda
# DISPLAY: RESMAX
# DISPLAY: RESMIN
# DISPLAY: ALAT
# DISPLAY: radlim
# DISPLAY: treatspotscan
# DISPLAY: ctfplotresmax
# DISPLAY: tempkeep
# DISPLAY: ISTEP
# DISPLAY: ISTEP_h
# DISPLAY: IMAXCOR
# DISPLAY: comment
# DISPLAY: RMAG
# DISPLAY: LCOLOR
# DISPLAY: createmaskinfo
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set PHASEORI_done = ""
set boxa1 = ""
set boxa2 = ""
set lattice = ""
set quadrada = ""
set imagename = ""
set nonmaskimagename = ""
set defocus = ""
set imagenumber = ""
set realcell = ""
set maska = ""
set imagesidelength = ""
set magnification = ""
set stepdigitizer = ""
set tempkeep = ""
set RESMIN = ""
set RESMAX = ""
set ALAT = ""
set quadpreda = ""
set radlim = ""
set realang = ""
set treatspotscan = ""
set facthresha = ""
set phacon = ""
set CS = ""
set KV = ""
set TLTAXIS = ""
set TLTANG = ""
set TLTAXA = ""
set TAXA = ""
set TANGL = ""
set ctfplotresmax = ""
set ISTEP = ""
set ISTEP_h = ""
set IMAXCOR = ""
set det_tilt = ""
set RMAG = ""
set LCOLOR = ""
set createmaskinfo = ""
set ctfcor_imode = ""
#
#$end_vars
#
set scriptname = 2dx_unbend1
\rm -f LOGS/${scriptname}.results
#
echo "<<@evaluate>>"
#
source ${proc_2dx}/initialize
#
if ( ${ctfcor_imode}x == 9x ) then
  set iname = image_ctfcor_multiplied
else
  set iname = image_ctfcor
endif
#
echo "<<@progress: 1>>"
#
set date = `date`
echo date = ${date}
#
if ( ! -e ${iname}.mrc ) then
  ${proc_2dx}/protest "Image ${iname}.mrc does not exist."
endif
#
if ( x${TLTAXIS} == "x-" || x${TLTAXIS} == "x--" ) then
  set TLTAXIS = "0.0"
endif
if ( x${TLTANG} == "x-" || x${TLTANG} == "x--" ) then
  set TLTANG = "0.0"
endif
if ( x${TLTAXA} == "x-" || x${TLTAXA} == "x--" ) then
  set TLTAXA = "0.0"
endif
if ( x${TAXA} == "x-" || x${TAXA} == "x--" ) then
  set TAXA = "0.0"
endif
if ( x${TANGL} == "x-" || x${TANGL} == "x--" ) then
  set TANGL = "0.0"
endif
#
#################################################################################
${proc_2dx}/linblock "Verifying some parameters"
#################################################################################
#
set tmp1 = `echo ${boxa1} | awk '{s = int( $1 ) } END { print s }'`
if ( ${tmp1} == ${boxa1} ) then
  echo boxa1 = ${boxa1}
else
  set boxa1 = ${tmp1}
  echo boxa1 = ${boxa1}
  echo "set boxa1 = ${boxa1}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Warning: boxa1 needs to be an integer number. Now corrected." 
  echo "#WARNING: Warning: boxa1 needs to be an integer number. Now corrected." >> LOGS/${scriptname}.results
endif

if ( ${boxa1} == "0" ) then
  set boxa1 = `echo ${imagesidelength} | awk '{ s = int( $1 / 10 ) } END { print s }'`
  echo "set boxa1 = ${boxa1}"  >> LOGS/${scriptname}.results
  echo ":Initializing boxa1 as ${boxa1}"
endif

set tmp1 = `echo ${boxa2} | awk '{s = int( $1 ) } END { print s }'`
if ( ${tmp1} == ${boxa2} ) then
  echo boxa2 = ${boxa2}
else
  set boxa2 = ${tmp1}
  echo boxa2 = ${boxa2}
  echo "set boxa2 = ${boxa2}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Warning: boxa2 needs to be an integer number. Now corrected."
  echo "#WARNING: Warning: boxa2 needs to be an integer number. Now corrected." >> LOGS/${scriptname}.results
endif

if ( ${boxa2} == "0" ) then
  set boxa2 = `echo ${imagesidelength} | awk '{ s = int( $1 / 11 ) } END { print s }'`
  echo "set boxa2 = ${boxa2}"  >> LOGS/${scriptname}.results
  echo ":Initializing boxa2 as ${boxa2}"
endif


set tmp1 = `echo ${boxa1} ${boxa2} | awk '{if ( $1 > $2 ) { s = 1 } else { s = 0 } } END { print s }'`
if ( ${tmp1} == 1 ) then
  set tmp1 = ${boxa2}
  set boxa2 = ${boxa1}
  set boxa1 = ${tmp1}
  echo boxa1 = ${boxa1}
  echo boxa2 = ${boxa2}
  echo "set boxa1 = ${boxa1}" >> LOGS/${scriptname}.results
  echo "set boxa2 = ${boxa2}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Warning: boxa2 should be larger than boxa1. Now exchanged."
  echo "#WARNING: Warning: boxa2 should be larger than boxa1. Now exchanged." >> LOGS/${scriptname}.results
endif
#

set rmax = 11000
echo rmax = ${rmax}
#
set imagecenterx = `echo ${imagesidelength} | awk '{ s = int( $1 / 2 ) } END { print s }'`
set imagecentery = ${imagecenterx}
#
#
set rtempx1 = ${imagecenterx}
set rtempy1 = ${imagecentery}
set rtempx2 = ${imagecenterx}
set rtempy2 = ${imagecentery}
@ rtempx1 -= 200
@ rtempx2 += 199
@ rtempy1 -= 200
@ rtempy2 += 199
# this gives a box at the reference locations with a diameter of 400 pixels.
set boxlabel = ${rtempx1},${rtempx2},${rtempy1},${rtempy2}
echo boxlabel = ${boxlabel}
#
set rtempx1 = ${imagecenterx}
set rtempy1 = ${imagecentery}
set rtempx2 = ${imagecenterx}
set rtempy2 = ${imagecentery}
@ rtempx1 -= 13
@ rtempx2 += 12
@ rtempy1 -= 13
@ rtempy2 += 12
# this gives a box at the reference location with a diameter of 26 pixels.
set patlabel = ${rtempx1},${rtempx2},${rtempy1},${rtempy2}
echo patlabel = ${patlabel}
#
# set RESPLOTMAX = 0.3
set RESPLOTMAX = `echo ${ctfplotresmax} | awk '{ if ( $1 > 0.1 ) { s = 1.0 / $1 } else { s = 0.3 } } END { print s }'`
echo RESPLOTMAX = ${RESPLOTMAX}
echo " 0.333 corresponds to 3.0 Angstroem for the border of the plot."
#
set quadradax = `echo ${quadrada} | sed 's/,/ /g' | awk '{ s = int( $1 ) } END { print s }'`
set quadraday = `echo ${quadrada} | sed 's/,/ /g' | awk '{ s = int( $2 ) } END { print s }'`
set tmp = `echo ${quadradax} ${quadraday} | awk '{if ( $2 == 0 ) { s = $1 } else { s = $2 }} END { print s }'`
if ( ${tmp} != ${quadraday} ) then
  set quadraday = ${tmp}
  set quadrada = `echo ${quadradax} ${quadraday} | sed 's/ /,/g'`
  echo "set quadrada = ${quadrada}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Warning: correcting quadrada to ${quadrada}."
  echo "#WARNING: Warning: correcting quadrada to ${quadrada}." >> LOGS/${scriptname}.results
endif
#
set u1 = `echo ${lattice} | cut -d\, -f1`
set u2 = `echo ${lattice} | cut -d\, -f2`
set v1 = `echo ${lattice} | cut -d\, -f3`
set v2 = `echo ${lattice} | cut -d\, -f4`
set ulen = `echo ${u1} ${u2} | awk '{ s = sqrt ( $1 * $1 + $2 * $2 ) } END { print s }'`
set vlen = `echo ${v1} ${v2} | awk '{ s = sqrt ( $1 * $1 + $2 * $2 ) } END { print s }'`
set smaller_vector = `echo ${ulen} ${vlen} | awk '{ if ( $1 < $2 ) { s = $1 } else { s = $2 }} END { print s }'`
#
set latlenok = `echo ${ulen} ${vlen} | awk '{ if ( $1 + $2 < 0.1 ) { s = 0 } else { s = 1 }} END { print s }'`
if ( ${latlenok} == '0' ) then
  ${proc_2dx}/linblock "ERROR: Lattice is ${lattice}"
  ${proc_2dx}/protest "Determine lattice first. Aborting."
endif
#
set newmaska = `echo ${smaller_vector} ${maska} | awk '{ if ( $1 / 2.1 < $2 ) { s = int( $1 / 2.1 ) } else { s = $2 }} END { print s }'`
if ( ${maska} != ${newmaska} ) then
  ${proc_2dx}/linblock "correcting maska from ${maska} to ${newmaska}"
  set maska = ${newmaska}
  echo "set maska = ${newmaska}" >> LOGS/${scriptname}.results
endif
${proc_2dx}/linblock "Using maska=${maska} for lattice of length u=${ulen} and v=${vlen}"
#
# Limit RADLIM parameters to at least 20
set radlimx = `echo ${radlim} | cut -d\, -f1`
set radlimy = `echo ${radlim} | cut -d\, -f2`
set radlima = `echo ${radlim} | cut -d\, -f3`
set radlimtmpx = `echo ${radlimx} | awk '{if ( $1 < 20 ) { s = 20 } else { s = $1 }} END { print s }'`
set radlimtmpy = `echo ${radlimy} | awk '{if ( $1 < 20 ) { s = 20 } else { s = $1 }} END { print s }'`
set radlimnewx = `echo ${radlimtmpx} | awk '{if ( $1 > 49 ) { s = 49 } else { s = $1 }} END { print s }'`
set radlimnewy = `echo ${radlimtmpy} | awk '{if ( $1 > 49 ) { s = 49 } else { s = $1 }} END { print s }'`
set radlimnew = `echo ${radlimnewx},${radlimnewy},${radlima}`
if ( ${radlim} != ${radlimnew} ) then
  set radlim = ${radlimnew}
  echo "set radlim = ${radlim}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "WARNING: correcting RADLIM to ${radlim}."
endif
#
echo "# IMAGE-IMPORTANT: "FFTIR/${iname}_fft.mrc "<FFT of Image>" >> LOGS/${scriptname}.results
echo "# IMAGE-IMPORTANT: "FFTIR/${iname}_red_fft.mrc "<FFT of Downsampled Image>" >> LOGS/${scriptname}.results
#
set final_round = "y"
#
#################################################################################
#
set outputfile = 2dx_unbend1.tmp
\rm -f ${outputfile}
#
${bin_2dx}/2dx_calcmag.exe << eot
${realcell}
${realang}
${TLTAXIS}
${TLTANG}
${lattice}
${imagesidelength}
${magnification}
${stepdigitizer}
${outputfile}
eot
#
if ( ! -e ${outputfile} ) then
  ${proc_2dx}/protest "ERROR in 2dx_calcmag.exe"
endif
#
set theormag = `cat ${outputfile} | head -n 1`
\rm -f ${outputfile}
#
${proc_2dx}/linblock "Theoretical magnification is ${theormag}, given magnification is ${magnification}"
#
echo "set CALCULATEDMAG = ${theormag}" >> LOGS/2dx_unbend1.results
#
if ( ! -e ${nonmaskimagename}.spt ) then
  echo ":: "
  ${proc_2dx}/linblock "Spotlist ${nonmaskimagename}.spt missing."
  ${proc_2dx}/protest "Aborting."
endif
#
if ( ! -e FFTIR/${iname}_fft.mrc ) then
  ${proc_2dx}/protest "ERROR: FFTIR/${iname}_fft.mrc not existing."
endif
#
echo "<<@progress: 5>>"
#
#################################################################################
${proc_2dx}/linblock "Sourcing 2dx_unbend1_sub.com"
#################################################################################
#
# this is only done for 2dx_unbend1_sub.com, so that it can stay as before:
set imagename = ${iname}
#
source ${proc_2dx}/2dx_unbend1_sub.com
#
#################################################################################
${proc_2dx}/linblock "Back in 2dx_unbend1.com"
#################################################################################
#
if ( ${ctfcor_imode}x == "x" || ${ctfcor_imode}x == "0x" ) then
  #
  #############################################################################
  #                                                                           #
  ${proc_2dx}/linblock "CTFAPPLY - to produce a first CTF overview"
  #                                                                           #
  #                     ${iname}_limit.aph  =>  ${iname}-limit-ctf.ps #
  #                                                                           #
  #############################################################################
  #
  \rm -f CTFPLOT.PS
  setenv IN  APH/${iname}_fou_unbend1.aph
  setenv OUT APH/${iname}_fou_unbend1_ctf.aph
  ${bin_2dx}/2dx_ctfapplyk.exe << eot
${lattice},${imagesidelength},${stepdigitizer},${magnification} ! AX,AY,BX,BY,ISIZE,DSTEP,XMAG
${defocus},${CS},${KV},${RESPLOTMAX} ! DFMID1,DFMID2,ANGAST,CS,KV,RESPLOTMAX
${imagenumber} ${iname}.u1, No Lim, Nyq=${ctfplotresmax}, ${date}
${phacon}
${RESMIN},0.0
${ctfcor_imode}  !  Mode of CTF correction
eot
  if ( ! -e CTFPLOT.PS ) then
    ${proc_2dx}/protest "ERROR in ctfapplyk.exe"
  else
    \mv -f CTFPLOT.PS PS/${iname}_fou_unbend1_ctf.ps
    echo "# IMAGE-IMPORTANT: PS/${iname}_fou_unbend1_ctf.ps <PS: IQ Plot>"  >> LOGS/${scriptname}.results
  endif
  #
  echo "<<@progress: 90>"
  #
#
endif
#
\rm -f fort.3
#
echo "<<@progress: 100>"
#
#############################################################################
${proc_2dx}/linblock "${scriptname} - normal end."
#############################################################################
#
#
